﻿using Alea;
using Alea.CSharp;
using Alea.CudaToolkit;
using Alea.Parallel;
using GSI.Core.Data;
using ManagedCuda;
using ManagedCuda.CudaBlas;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
using OpenCL.Net;
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace GSI.Core
{
	abstract class KriGing
	{
		public KriGing() { }

		protected class RegressionPoint
		{
			public double distance;
			public double Semivariance;

			public RegressionPoint(double distance, double semivariogram)
			{
				this.distance = distance;
				this.Semivariance = semivariogram;
			}
			public RegressionPoint()
			{ }
		}

		protected List<GPoint> pointList;

		protected double formula_a;
		protected double formula_r;
		protected int size;

		/// <summary>
		/// 构造函数 根据传入的采样点计算拟合曲线
		/// </summary>
		/// <param name="pointList"></param>
		protected KriGing(List<GPoint> pointList)
		{
			this.pointList = pointList;
			size = pointList.Count;

			#region 计算 formula_c formula_r 拟合曲线"f(x)=a*(1-e^(-x/r))"中的a，r   
			List<RegressionPoint> initialPoints = new List<RegressionPoint>();

			for (int m = 0; m < pointList.Count; m++)
				for (int n = 0; n < pointList.Count; n++)
					if (m != n)
						initialPoints.Add(GetRegressionPoint(pointList[m], pointList[n]));

			double maxDistance = initialPoints.Max(e => e.distance);
			double minDistance = initialPoints.Min(e => e.distance);

			List<RegressionPoint> showPoints = new List<RegressionPoint>();
			List<RegressionPoint> calPoints = new List<RegressionPoint>();

			for (int n = (int)minDistance; n <= maxDistance; n++)
			{
				//此处tempPointsList为0
				var tempPointsList = initialPoints.Where(e => (int)e.distance == n);

				if (tempPointsList.Count() > 0)
				{
					showPoints.Add(new RegressionPoint {
						distance = n,
						Semivariance = tempPointsList.Average(e => e.Semivariance) // 平均值
					});
				}
			}

			for (int n = 0; n < showPoints.Count / 10; n++)
			{
				var tmp = showPoints.GetRange(n * 10, 10);

				calPoints.Add(new RegressionPoint {
					distance = tmp.Average(p => p.distance),
					Semivariance = tmp.Average(p => p.Semivariance)
				});
			}

			//高斯牛顿法
			pgPowerModel model = new pgPowerModel();
			pgGaussNewtonSolver solver = new pgGaussNewtonSolver(0.001, 0.001, 10000, new DenseVector(new[] { 50.0, 1.5 }));
			List<Vector<double>> solverIterations = new List<Vector<double>>();

			double[] x = new double[calPoints.Count];
			double[] y = new double[calPoints.Count];
			for (int n = 0; n < calPoints.Count; n++)
			{
				//缩小数值，方便计算，之后会还原
				x[n] = calPoints[n].distance / 5;
				y[n] = calPoints[n].Semivariance / 1000;
			}
			calPoints.ForEach(e => { e.distance /= 5; e.Semivariance /= 1000; });

			solver.Estimate(model, calPoints.Count, new DenseVector(x), new DenseVector(y), ref solverIterations);

			formula_a = solverIterations.Last()[0] * 5;
			formula_r = solverIterations.Last()[1] * 1000;

			#endregion
		}

		protected void InitKriGing(List<GPoint> pointList)
		{
			this.pointList = pointList;
			size = pointList.Count;

			#region 计算 formula_c formula_r 拟合曲线"f(x)=a*(1-e^(-x/r))"中的a，r   
			List<RegressionPoint> initialPoints = new List<RegressionPoint>();

			for (int m = 0; m < pointList.Count; m++)
				for (int n = 0; n < pointList.Count; n++)
					if (m != n)
						initialPoints.Add(GetRegressionPoint(pointList[m], pointList[n]));

			double maxDistance = initialPoints.Max(e => e.distance);
			double minDistance = initialPoints.Min(e => e.distance);

			List<RegressionPoint> showPoints = new List<RegressionPoint>();
			List<RegressionPoint> calPoints = new List<RegressionPoint>();

			for (int n = (int)minDistance; n <= maxDistance; n++)
			{
				//此处tempPointsList为0
				var tempPointsList = initialPoints.Where(e => (int)e.distance == n);

				if (tempPointsList.Count() > 0)
				{
					showPoints.Add(new RegressionPoint {
						distance = n,
						Semivariance = tempPointsList.Average(e => e.Semivariance) // 平均值
					});
				}
			}

			for (int n = 0; n < showPoints.Count / 10; n++)
			{
				var tmp = showPoints.GetRange(n * 10, 10);

				calPoints.Add(new RegressionPoint {
					distance = tmp.Average(p => p.distance),
					Semivariance = tmp.Average(p => p.Semivariance)
				});
			}

			//高斯牛顿法
			pgPowerModel model = new pgPowerModel();
			pgGaussNewtonSolver solver = new pgGaussNewtonSolver(0.001, 0.001, 10000, new DenseVector(new[] { 50.0, 1.5 }));
			List<Vector<double>> solverIterations = new List<Vector<double>>();

			double[] x = new double[calPoints.Count];
			double[] y = new double[calPoints.Count];
			for (int n = 0; n < calPoints.Count; n++)
			{
				//缩小数值，方便计算，之后会还原
				x[n] = calPoints[n].distance / 5;
				y[n] = calPoints[n].Semivariance / 1000;
			}
			calPoints.ForEach(e => { e.distance /= 5; e.Semivariance /= 1000; });

			solver.Estimate(model, calPoints.Count, new DenseVector(x), new DenseVector(y), ref solverIterations);

			formula_a = solverIterations.Last()[0] * 5;
			formula_r = solverIterations.Last()[1] * 1000;

			#endregion
		}
		/// <summary>
		/// 克里金插值
		/// </summary>
		/// <param name="meshPointList">网格点列表</param>
		/// <returns>存有插值数据的网格点列表</returns>
		abstract public string GetInterpolationData(List<GPoint> meshPointList);

		public Info info = new Info();

		/// <summary>
		/// 根据拟合曲线计算半方差
		/// </summary>
		protected double GetRij(double x1, double y1, double x2, double y2)
		{
			double distance = System.Math.Sqrt(System.Math.Pow(x1 - x2, 2) + System.Math.Pow(y1 - y2, 2));
			return distance == 0 ? formula_a : formula_a * System.Math.Exp(-distance / formula_r);
		}
		/// <summary>
		/// 根据拟合曲线计算半方差
		/// </summary>
		protected double GetRij(GPoint p1, GPoint p0)
		{
			double distance = GetDistance(p1,p0);
			return distance == 0 ? formula_a : formula_a * System.Math.Exp(-distance / formula_r);
		}
		/// <summary>
		/// 根据 1/2 （z1 - z0）^2 计算两点半方差
		/// </summary>
		protected double GetSemivariance(GPoint p1, GPoint p0) => 0.5 * System.Math.Pow((p1.Value - p0.Value), 2);

		/// <summary>
		/// 根据两点计算距离
		/// </summary>
		protected double GetDistance(double x1, double y1, double x0, double y0) => System.Math.Sqrt(System.Math.Pow(x1 - x0, 2) + System.Math.Pow(y1 - y0, 2));
		/// <summary>
		/// 根据两点计算距离
		/// </summary>
		protected double GetDistance(GPoint p1, GPoint p0) => GetDistance(p1.X, p1.Y, p0.X, p0.Y);
		/// <summary>
		/// 根据两点生成回归点
		/// </summary>
		protected RegressionPoint GetRegressionPoint(GPoint p1, GPoint p0) => new RegressionPoint {
			distance = GetDistance(p1, p0),
			Semivariance = GetSemivariance(p1, p0)
		};
	}

	/// <summary>
	/// 使用 串行 计算
	/// </summary>
	class CPUSerial : KriGing
	{
		public CPUSerial(List<GPoint> pointList) : base(pointList)
		{
		}

		public override string GetInterpolationData(List<GPoint> meshPointList)
		{
			string output = "**************************\n";
			Stopwatch sw = new Stopwatch();
			sw.Start();

			//计算协方差 取逆
			var k = new double[size, size];
			for (int m = 0; m < size; m++)
				for (int n = 0; n < size; n++)
				{
					k[m, n] = GetRij(pointList[m], pointList[n]);
				}

			var kInverse = CreateMatrix.DenseOfArray<double>(k).LU().Inverse().ToArray();
			sw.Stop();
			var t0 = (sw.ElapsedMilliseconds);
			output += $"计算协方差时间{t0}ms\n";
			sw.Start();
			// 未知点的半方差
			double[,] Dm = new double[pointList.Count, meshPointList.Count];

			for (int m = 0; m < meshPointList.Count; m++)
			{
				for (int p = 0; p < size; p++)
					Dm[p, m] = GetRij(pointList[p], meshPointList[m]);
			}
			sw.Stop();
			var t1 = (sw.ElapsedMilliseconds);
			output += $"计算未知点的半方差时间{t1 - t0}ms\n";
			sw.Start();
			// 串行计算
			double[,] namuta = new double[size, meshPointList.Count];
			CPUSerial.RunCpu(kInverse, Dm, namuta);
			sw.Stop();
			var t2 = (sw.ElapsedMilliseconds);
			output += $"计算矩阵相乘时间{t2 - t1}ms\n";
			sw.Start();
			// 网格点估值
			for (int m = 0; m < meshPointList.Count; m++)
			{
				for (int q = 0; q < size; q++)
					meshPointList[m].Value += namuta[q, m] * pointList[q].Value;
			}
			sw.Stop();
			var t3 = (sw.ElapsedMilliseconds);
			output += $"计算累加估值时间{t3 - t2}ms\n克里金串行计算总时间{t3}ms\n***************************\n";
			info.Time_Km = t0;
			info.Time_Dm = t1 - t0;
			info.Time_MM = t2 - t1;
			info.Time_Re = t3 - t2;
			info.Time_All = t3;
			return output;
		}

		private static void RunCpu(double[,] a, double[,] b, double[,] c)
		{
			for (var row = 0; row < c.GetLength(0); ++row)
			{
				for (var col = 0; col < c.GetLength(1); ++col)
				{
					var sum = 0.0;
					for (var k = 0; k < a.GetLength(1); ++k)
					{
						sum += a[row, k] * b[k, col];
					}
					c[row, col] = sum;
				}
			}
		}
	}

	/// <summary>
	/// 使用 多核 实现克里金算法
	/// </summary>
	class CPUParallel : KriGing
	{
		public CPUParallel() { }
		public CPUParallel(List<GPoint> pointList) : base(pointList) { }
		public override string GetInterpolationData(List<GPoint> meshPointList)
		{
			string output = "**************************\n";
			Stopwatch sw = new Stopwatch();
			sw.Start();

			//计算协方差 取逆
			var k = new double[size, size];
			Parallel.For(0, size, (m) =>
			{
				for (int n = 0; n < size; n++)
				{
					k[m, n] = GetRij(pointList[m], pointList[n]);
				}
			});

			var klu = CreateMatrix.DenseOfArray<double>(k).LU().Inverse();
			sw.Stop();
			var t0 = (sw.ElapsedMilliseconds);
			output += $"计算协方差时间{t0}ms\n";
			sw.Start();


			// 修改处 多线程使用[]
			Parallel.For(0, meshPointList.Count, (m) =>
			{
				var Dm = new DenseVector(size);
				for (int p = 0; p < size; p++)
					Dm[p] = GetRij(pointList[p], meshPointList[m]);
				var namuta = klu.LeftMultiply(Dm);
				for (int q = 0; q < size; q++)
				{
					meshPointList[m].Value += namuta[q] * pointList[q].Value;
				}
			});
			sw.Stop();
			var t3 = (sw.ElapsedMilliseconds);
			output += $"克里金多核多线程计算总时间{t3}ms\n***************************\n";
			info.Time_Km = t0;
			//info.Time_Dm = t1 - t0;
			//info.Time_MM = t2 - t1;
			//info.Time_Re = t3 - t2;
			info.Time_All = t3;
			return output;
		}
	}

	/// <summary>
	/// 使用 OpencL 实现克里金算法
	/// </summary>
	class GPUOpenCL : KriGing
	{
		public GPUOpenCL(List<GPoint> pointList) : base(pointList) { }
		public override string GetInterpolationData(List<GPoint> meshPointList)
		{

			string output = "**************************\n";
			Stopwatch sw = new Stopwatch();

			//计算协方差 取逆
			sw.Start();
			var k = new double[size, size];
			Parallel.For(0, size, (m) =>
			{
				for (int n = 0; n < size; n++)
				{
					k[m, n] = GetRij(pointList[m], pointList[n]);
				}
			});

			var kInverse = CreateMatrix.DenseOfArray<double>(k).LU().Inverse();
			sw.Stop();

			var t0 = (sw.ElapsedMilliseconds);
			output += $"计算协方差时间{t0}ms\n";

			// 计算半方差
			sw.Start();
			double[] Dm = new double[size*meshPointList.Count];
			Parallel.For(0, meshPointList.Count, (m) =>
			{
				for (int p = 0; p < size; p++)
					Dm[p * meshPointList.Count + m] = GetRij(pointList[p], meshPointList[m]);
			});
			sw.Stop();

			var t1 = (sw.ElapsedMilliseconds);
			output += $"计算未知点的半方差时间{t1 - t0}ms\n";

			// 矩阵
			sw.Start();
			float[] namuta = OpenclTools.GPUMartixMultiply(pointList.Count, pointList.Count, meshPointList.Count, kInverse, Dm);
			sw.Stop();

			var t2 = (sw.ElapsedMilliseconds);
			output += $"计算矩阵相乘时间{t2 - t1}ms\n";

			// 插值
			sw.Start();
			Parallel.For(0, meshPointList.Count, (m) =>
			{
				for (int q = 0; q < size; q++)
					meshPointList[m].Value += namuta[q + m * size] * pointList[q].Value;
			});
			sw.Stop();
			var t3 = (sw.ElapsedMilliseconds);
			output += $"计算累加估值时间{t3 - t2}ms\n克里金OpenCL计算总时间{t3}ms\n***************************\n";
			info.Time_Km = t0;
			info.Time_Dm = t1 - t0;
			info.Time_MM = t2 - t1;
			info.Time_Re = t3 - t2;
			info.Time_All = t3;
			return output;
		}
	}

	class OpenclTools
	{
		private static CommandQueue commandQueue;
		private static OpenCL.Net.Context context;
		private static ErrorCode error;
		private static OpenCL.Net.Kernel kernel;
		private static OpenCL.Net.Program program;
		private static bool init = false;
		private static string source = @"
            //Matrix Multiply
            __kernel void mul_matrix(
                const int iN,
                const int iM,
                const int iP,
                __global const float* A, 
                __global const float* B, 
                __global float* C)
            {
                int i = get_global_id(0);
                int k,j;
                float tmp;
                if ((i < iN)) {
                    for(j = 0; j < iM; j++)
                    {
                        tmp = 0.0;
                        for (k = 0; k < iP; k++)
                            tmp += A[i*iP + k] * B[k*iM + j];
                        C[i*iM + j] = tmp;
                    }
                }
            }";
		private static ulong kernel_start_to_end_time;

		static OpenclTools()
		{
			if (init)
				return;
			// 1. Allocates memory for host buffers and initializes them.
			// 已经完成

			// 2. Gets platform and device information. 
			//    This is discussed in detail in Chapter 2, OpenCL Architecture.
			Platform[] platforms = Cl.GetPlatformIDs(out error);
			// 3. Sets up the platform.
			Platform platform = platforms[0];
			// 4. Gets the devices list and chooses the type of device you want to run on.
			OpenCL.Net.Device[] devices = Cl.GetDeviceIDs(platforms[0], DeviceType.Gpu, out error);
			// Log.WriteLine("使用设备信息：");
			foreach (var item in devices)
			{
				var str = Cl.GetDeviceInfo(item, DeviceInfo.Name, out error);
				// Log.WriteLine(str.ToString());
			}
			OpenCL.Net.Device device = devices[0];
			// 5. Creates an OpenCL context for the device.
			context = Cl.CreateContext(null, 1, devices, null, IntPtr.Zero, out error);
			// 6. Creates a command queue.
			commandQueue = Cl.CreateCommandQueue(context, device, CommandQueueProperties.ProfilingEnable, out error);
			// 9. Creates a program from the kernel source.
			program = Cl.CreateProgramWithSource(context, 1, new[] { source }, null, out error);
			// 10. Builds the program and creates the OpenCL kernel.
			Cl.BuildProgram(program, 0, null, string.Empty, null, IntPtr.Zero);
			// check build
			if (Cl.GetProgramBuildInfo(program, device, ProgramBuildInfo.Status, out error).CastTo<BuildStatus>() != BuildStatus.Success)
			{
				if (error != ErrorCode.Success)
					Console.WriteLine("ERROR: " + "Cl.GetProgramBuildInfo" + " (" + error.ToString() + ")");
				Console.WriteLine("Cl.GetProgramBuildInfo != Success");
				Console.WriteLine(Cl.GetProgramBuildInfo(program, device, ProgramBuildInfo.Log, out error));
			}
			kernel = Cl.CreateKernel(program, "mul_matrix", out error);
			init = true;
		}

		public static float[] GPUMartixMultiply(int N, int P, int M, Matrix<double> kn, double[] dm)
		{

			#region 数据初始化
			float[] A;
			float[] B;
			float[] C;

			A = new float[N * P];
			B = new float[P * M];
			C = new float[N * M];

			for (int i = 0; i < N; i++)
			{
				for (int j = 0; j < P; j++)
				{
					A[i * P + j] = (float)kn.At(i, j);
				}
			}
			for (int i = 0; i < P * M; i++)
			{
				B[i] = (float)dm[i];
			}
			#endregion

			#region GPU 并行执行

			// 7. Creates memory buffers on the device for each vector.
			Mem mem_A = (Mem)Cl.CreateBuffer(context, MemFlags.ReadOnly, sizeof(float) * A.Count(), out error);
			Mem mem_B = (Mem)Cl.CreateBuffer(context, MemFlags.ReadOnly, sizeof(float) * B.Count(), out error);
			Mem mem_C = (Mem)Cl.CreateBuffer(context, MemFlags.ReadWrite, sizeof(float) * C.Count(), out error);
			// 8. Copies the Buffer A and B to the device.
			OpenCL.Net.Event event0;
			OpenCL.Net.Event event_kernel;
			Cl.EnqueueWriteBuffer(commandQueue, mem_A, Bool.True, IntPtr.Zero, new IntPtr(sizeof(float) * A.Count()), A, 0, null, out event0);
			Cl.EnqueueWriteBuffer(commandQueue, mem_B, Bool.True, IntPtr.Zero, new IntPtr(sizeof(float) * B.Count()), B, 0, null, out event0);
			Cl.EnqueueWriteBuffer(commandQueue, mem_C, Bool.True, IntPtr.Zero, new IntPtr(sizeof(float) * C.Count()), C, 0, null, out event0);

			// 11. Sets the arguments of the kernel.
			error = Cl.SetKernelArg(kernel, 0, N);
			error = Cl.SetKernelArg(kernel, 1, M);
			error = Cl.SetKernelArg(kernel, 2, P);
			error = Cl.SetKernelArg(kernel, 3, mem_A);
			error = Cl.SetKernelArg(kernel, 4, mem_B);
			error = Cl.SetKernelArg(kernel, 5, mem_C);

			// 12. Executes the OpenCL kernel on the device.
			error = Cl.EnqueueNDRangeKernel(commandQueue, kernel, 1, null, new IntPtr[] { new IntPtr(N) }, null, 0, null, out event_kernel);
			// 13. Reads back the memory from the device to the host buffer. 
			//  This step is optional, you may want to keep the data resident 
			//      in the device for further processing.
			error = Cl.EnqueueReadBuffer(commandQueue, mem_C, Bool.True, 0, C.Count(), C, 0, null, out event0);
			// 14. Cleans up and waits for all the commands to complete.
			Cl.Finish(commandQueue);

			Cl.WaitForEvents(1, new OpenCL.Net.Event[] { event_kernel });
			InfoBuffer vara, varb;
			vara = Cl.GetEventProfilingInfo(event_kernel, ProfilingInfo.Start, out error); //OpenCL.Net.ErrorCode.InvalidEvent
			varb = Cl.GetEventProfilingInfo(event_kernel, ProfilingInfo.End, out error);
			kernel_start_to_end_time = (varb.CastTo<ulong>() - vara.CastTo<ulong>());

			// 15. Finally releases all OpenCL allocated objects and host buffers.
			Cl.ReleaseMemObject(mem_A);
			Cl.ReleaseMemObject(mem_B);
			Cl.ReleaseMemObject(mem_C);
			#endregion

			return C;
		}

		public static float[] CublasMartixMultiply(int N, int P, int M, Matrix<double> kn, double[] dm)
		{
			return null;
		}
	}

    /// <summary>
    /// 使用 ManagedCuda- Cublas 实现克里金算法
    /// </summary>
    class GPUMCublas : KriGing
    {
        public GPUMCublas(List<GPoint> pointList) : base(pointList) { }
        public override string GetInterpolationData(List<GPoint> meshPointList)
        {
            string output = "**************************\n";
            Stopwatch sw = new Stopwatch();
            sw.Start();

            //计算协方差 取逆
            var k = new double[size, size];
            double[] Km = new double[size * size];
            Parallel.For(0, size, (m) =>
            {
                for (int n = 0; n < size; n++)
                {
                    k[m, n] = GetRij(pointList[m], pointList[n]);
                }
            });

            var kInverse = CreateMatrix.DenseOfArray<double>(k).LU().Inverse().ToArray();
            // 降为1为数组
            Parallel.For(0, size, (i) =>
            {
                for (int j = 0; j < size; j++)
                {
                    Km[i * size + j] = kInverse[i, j];
                }
            });

            sw.Stop();
            var t0 = (sw.ElapsedMilliseconds);
            output += $"计算协方差时间{t0}ms\n";
            sw.Start();

            // 网格点插值
            double[] Dm = new double[size * meshPointList.Count];
            Parallel.For(0, meshPointList.Count, (m) =>
            {
                for (int p = 0; p < size; p++)
                    Dm[p * meshPointList.Count + m] = GetRij(pointList[p], meshPointList[m]);
            });

            sw.Stop();
            var t1 = (sw.ElapsedMilliseconds);
            output += $"计算未知点的半方差时间{t1 - t0}ms\n";
            sw.Start();
            double[] namuta = new double[size * meshPointList.Count];
			//TestMN();

			namuta = RunMCublas(size, size, meshPointList.Count(), Km, Dm);
            sw.Stop();
            var t2 = (sw.ElapsedMilliseconds);
            output += $"计算矩阵相乘时间{t2 - t1}ms\n";
            sw.Start();
            Parallel.For(0, meshPointList.Count, (m) =>
            {
                for (int q = 0; q < size; q++)
                    meshPointList[m].Value += namuta[q * meshPointList.Count + m] * pointList[q].Value;
            });
            sw.Stop();
            var t3 = (sw.ElapsedMilliseconds);
            output += $"计算累加估值时间{t3 - t2}ms\n克里金Cublas算总时间{t3}ms\n***************************\n";
			info.Time_Km = t0;
			info.Time_Dm = t1 - t0;
			info.Time_MM = t2 - t1;
			info.Time_Re = t3 - t2;
			info.Time_All = t3;
			return output;
        }

        public double[] RunMCublas(int M,int K,int N,double[] Km,double[] Dm)
        {
            int nA = M * K;
            int nB = K * N;
            int nC = M * N;

            double[] h_A = Km;
            double[] h_B = Dm;
            double[] h_C;

            CudaDeviceVariable<double> d_A;
            CudaDeviceVariable<double> d_B;
            CudaDeviceVariable<double> d_C;
            double alpha = 1.0;
            double beta = 0.0;
            CudaBlas handle;


            /* Initialize CUBLAS */
            Console.WriteLine("simpleCUBLAS test running.");
            handle = new CudaBlas();

            /* Allocate host memory for the matrices */
			

            /* Allocate device memory for the matrices */
            d_A = new CudaDeviceVariable<double>(nA);
            d_B = new CudaDeviceVariable<double>(nB);
            d_C = new CudaDeviceVariable<double>(nC);

            Stopwatch s = new Stopwatch();
            s.Start();
            /* Initialize the device matrices with the host matrices */
            d_A.CopyToDevice(h_A);
            d_B.CopyToDevice(h_B);
            //d_C.CopyToDevice(h_C);
            /* Performs operation using cublas */
            //handle.Gemm(Operation.NonTranspose, Operation.NonTranspose, M, N, K, alpha, d_A, K, d_B, N, beta, d_C, N); // wrong
            //handle.Gemm(Operation.Transpose, Operation.Transpose, M, N, K, alpha, d_A, K, d_B, N, beta, d_C, M);// right by col 返回A*B的转置
            handle.Gemm(Operation.NonTranspose, Operation.NonTranspose, N, M, K, alpha, d_B, N, d_A, K, beta, d_C, N); // right by row 
            s.Stop();

            /* Allocate host memory for reading back the result from device memory */
            h_C = d_C;

            /* Memory clean up */
            d_A.Dispose();
            d_B.Dispose();
            d_C.Dispose();

            /* Shutdown */
            handle.Dispose();
            Console.WriteLine("simpleCUBLAS test finished.");
			return h_C;
        }
		private void simple_sgemm(int mar_m, int mar_n, int mar_k, double alpha, double[] A, double[] B, double beta, double[] C)
		{
			int i;
			int j;
			int k;

			for (i = 0; i < mar_m; ++i)
			{
				for (j = 0; j < mar_n; ++j)
				{
					double prod = 0;

					for (k = 0; k < mar_k; ++k)
					{
						prod += A[i * mar_k + k] * B[k * mar_n + j];
					}

					C[i * mar_n + j] = alpha * prod + beta * C[i * mar_n + j];
				}
			}
		}
	}

	/// <summary>
	/// 使用 CUDA 实现克里金算法
	/// </summary>
	class GPUCUDA : KriGing
	{
		public GPUCUDA(List<GPoint> pointList) : base(pointList) { }
		public override string GetInterpolationData(List<GPoint> meshPointList)
		{
			string output = "**************************\n";
			
			Stopwatch sw = new Stopwatch();
			sw.Start();


			//计算协方差 取逆 转置
			var k = new double[size, size];
			double[,] Km = new double[size, size];
			Parallel.For(0, size, (m) =>
			{
				for (int n = 0; n < size; n++)
				{
					k[m, n] = GetRij(pointList[m], pointList[n]);
				}
			});

			var kInverse = CreateMatrix.DenseOfArray<double>(k).LU().Inverse().ToArray();
			sw.Stop();
			var t0 = (sw.ElapsedMilliseconds);
			output += $"计算协方差时间{t0}ms\n";
			sw.Start();

			// 网格点插值
			double[,] Dm = new double[pointList.Count, meshPointList.Count];
			Parallel.For(0, meshPointList.Count,
				(m) =>
				{
					for (int p = 0; p < size; p++)
						Dm[p, m] = GetRij(pointList[p], meshPointList[m]);
				});
			sw.Stop();
			var t1 = (sw.ElapsedMilliseconds);
			output += $"计算未知点的半方差时间{t1 - t0}ms\n";
			sw.Start();
			double[,] namuta = new double[size, meshPointList.Count];

			RunGpu(kInverse, Dm, namuta);
			sw.Stop();
			var t2 = (sw.ElapsedMilliseconds);
			output += $"计算矩阵相乘时间{t2 - t1}ms\n";
			sw.Start();
			Parallel.For(0, meshPointList.Count, (m) =>
			{
				for (int q = 0; q < size; q++)
					meshPointList[m].Value += namuta[q, m] * pointList[q].Value;
			});

			sw.Stop();
			var t3 = (sw.ElapsedMilliseconds);
			output += $"计算累加估值时间{t3 - t2}ms\n克里金CUDA计算总时间{t3}ms\n***************************\n";
			info.Time_Km = t0;
			info.Time_Dm = t1 - t0;
			info.Time_MM = t2 - t1;
			info.Time_Re = t3 - t2;
			info.Time_All = t3;
			return output; ;
		}

		#region CUDA 编程

		private const int BlockSize = 32;
		/// <summary>
		/// GPU计算 AB=C
		/// </summary>
		/// <param name="a"></param>
		/// <param name="b"></param>
		/// <param name="c"></param>
		[GpuManaged]
		private void RunGpu(double[,] a, double[,] b, double[,] c)
		{
			var lp = LaunchParam(a, b, c);
			Gpu.Default.Launch(Kernel, lp, a, b, c);
		}


		private double GetMatrixElement(double[,] matrix, int blockRow, int blockCol, int row, int col)
		{
			var globalRow = blockRow * BlockSize + row;
			var globalCol = blockCol * BlockSize + col;
			if (globalRow < matrix.GetLength(0) && globalCol < matrix.GetLength(1))
				return matrix[globalRow, globalCol];
			else
				return 0.0;
		}

		private void SetMatrixElement(double[,] matrix, int blockRow, int blockCol, int row, int col, double value)
		{
			var globalRow = blockRow * BlockSize + row;
			var globalCol = blockCol * BlockSize + col;
			if (globalRow < matrix.GetLength(0) && globalCol < matrix.GetLength(1))
				matrix[globalRow, globalCol] = value;
		}

		private int DivUp(int num, int den)
		{
			return (num + den - 1) / den;
		}

		private void Kernel(double[,] a, double[,] b, double[,] c)
		{
			var colsA = a.GetLength(1);
			var blockRow = blockIdx.x;
			var blockCol = blockIdx.y;

			var valueC = 0.0;

			var row = threadIdx.x;
			var col = threadIdx.y;

			for (var m = 0; m < DivUp(colsA, BlockSize); ++m)
			{
				var subA = __shared__.Array2D<double>(BlockSize, BlockSize);
				var subB = __shared__.Array2D<double>(BlockSize, BlockSize);

				subA[row, col] = GetMatrixElement(a, blockRow, m, row, col);
				subB[row, col] = GetMatrixElement(b, m, blockCol, row, col);
				DeviceFunction.SyncThreads();

				for (var e = 0; e < BlockSize; ++e)
				{
					valueC += subA[row, e] * subB[e, col];
				}
				DeviceFunction.SyncThreads();
			}

			SetMatrixElement(c, blockRow, blockCol, row, col, valueC);
		}


		private LaunchParam LaunchParam(double[,] a, double[,] b, double[,] c)
		{
			Check(a, b, c);
			var blockSize = new Alea.dim3(BlockSize, BlockSize);
			var gridSize = new Alea.dim3(DivUp(a.GetLength(0), BlockSize), DivUp(b.GetLength(1), BlockSize));
			return new LaunchParam(gridSize, blockSize);
		}

		private void Check(double[,] a, double[,] b, double[,] c)
		{
			if (a == null) throw new ArgumentNullException(nameof(a));
			if (b == null) throw new ArgumentNullException(nameof(b));
			if (c == null) throw new ArgumentNullException(nameof(c));
			Debug.Assert(a.GetLength(1) == b.GetLength(0));
			Debug.Assert(a.GetLength(0) == c.GetLength(0));
			Debug.Assert(b.GetLength(1) == c.GetLength(1));
		}
		#endregion
	}

}


